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The /-electron spectral function of the Falicov-Kimball model is calculated via a Keldysh-based 
many-body formalism originally developed by Brandt and Urbanek. We provide results for both the 
Bethe lattice and the hypercubic lattice at half filling. Since the numerical computations are quite 
sensitive to the discretization along the Kadanoff-Baym contour and to the maximum cutoff in time 
that is employed, we analyze the accuracy of the results using a variety of different moment sum-rules 
and spectral formulas. We find that the /-electron spectral function has interesting temperature 
dependence becoming a narrow single-peaked function for small U and developing a gap, with two 
broader peaks for large U. 
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I. INTRODUCTION 



Nonequilibrium many-body physics is becoming an in- 
creasingly important field because it allows strongly cor- 
related electrons to be examined in the presence of large 
external fields, which can drive them far from equilib- 
rium; as nanotechnology research grows, there are an in- 
creasing number of nanoscale solid-state devices which 
are exposed to extreme fields. The formalism to solve 
nonequilibrium problems was developed independently 
by KadanofF and Baym^ and by Keldysh^. In princi- 
ple, it allows these nonequilibrium problems to be solved 
exactly (including all nonlinear field effects), but usually 
the formalism is employed in a perturbative approach 
for the electron correlations. These nonequilibrium for- 
malisms can also be used to calculate equilibrium Green's 
functions, especially in cases where alternative analytic 
continuation methods are intractable. We examine the 
simplest such problem — the /-electron spectrum of the 
Falicov-Kimball model. This problem can be solved ex- 
actly with dynamical mean field theory (DMFT). 

The spinless Falicov-Kimball model? (FK) describes 
the dynamics of two types of electrons: conduction elec- 
trons (created or destroyed at site i by d'^ or d) and lo- 
calized electrons (created or destroyed at site i by /''' or 
/). The non- interacting conduction electrons can hop 
between nearest-neighbor sites on a D-dimensional lat- 
tice, with a hopping matrix —tij — —t*/2^/D; we choose 
a scaling of the hopping matrix that yields a nontrivial 
limit in infinite-dimensions"'. The f-electrons have a site 
energy Ef, and a chemical potential ^ is employed to ad- 
just the total number of electrons rid + Uf = '^tot- The 
d- and f- number operators at each site are Ud and rif. 
There is a finite Coulomb interaction U between d- and 
f-electrons that occupy the same lattice site, and so the 



Falicov-Kimball Hamiltonian for the lattice is 

ij i 

+ J^udldjju (1) 

i 

The FK lattice model [in Eq. lU] can be solved in infinite 
dimensions using the methods of Brandt and Mielsch^. 
We consider two kinds of lattices: (i) the hypercu- 
bic lattice with a Gaussian noninteracting density of 
states /9?tc(e) = exp[— e^/t*^]/(-y7rr), and (n) the infinite- 
coordination Bethe lattice with a semicircular noninter- 
acting density of states Pb{^) — V^t*'^ — e^/(27rt*^); we 
take t* as the unit of energy {t* = 1), and consider only 
the homogeneous phase, where all quantities are transla- 
tionally invariant. 

Despite the conservation of the local /-electron num- 
ber {[Hfk, fj , fi] = 0), the /-electrons have nontriv- 
ial dynamics as a function of T^'^, so we expect the 
/-electron spectral function to have an interesting evo- 
lution with the correlation strength U. When U = 0, 
the spectral function is a delta function, which broad- 
ens as U increases because the /-electron sees a fluctu- 
ating d-electron that hops onto and off of the local site 
as a function of time. When U is increased further, a 
metal-insulator transition takes place, and we expect the 
/-spectral function to develop a gap as T — > 0. Hence the 
/-electron spectral function should have rich behavior as 
a function of U and T. 

The organization of this paper is as follows. In Section 
II, we present the formalism, outlining in detail how the 
analytic continuation is carried out within this nonequi- 
librium approach; our notation is nonstandard because 
we use the usual Kadanoff-Baym contour, but we stick 
with the Green's function definitions for and of 
Brandt and Urbanek^. In Section III, we present our 
numerical results for both the Bethe and hypercubic lat- 
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tices, and we conclude in Section IV. A short communica- 
tion, showing the spectral function at the critical or near 
critical value of U for the metal-insulator transition has 
been completed^, and a longer contribution, detailing the 
computational algorithm, parallelization, and numerical 
accuracy, has also been completediii. 

II. FORMALISM 

The many body problem on an infinite-coordination 
lattice can be solved by a mean-field-like procedure, be- 
cause the self energy of the conduction electrons is local^. 
Hence, the local d-electron Green's function Gioc{z) on 
the lattice satisfies 

G'n^) - / ^ ^^'\ , de, (2) 

where z is a complex variable, and is the momentum- 
independent self energy. As noted by Brandt and 
MielscbSi, the lattice self energy coincides with the self 
energy of an atomic d-state coupled to an f-state with 
the same Coulomb interaction as on the lattice, and 
perturbed by an external time-dependent field, A(t, r'), 
which mimics the hopping on the lattice. For an appro- 
priate choice of the A-field, the functional dependence 
of Srf on Gd{z) and Gf{z), the atomic propagators for 
d- and f-states, is exactly the same as in the lattice case. 
The lattice problem is thus reduced to finding the atomic 
self-energy functional for the d-electrons, and then set- 
ting G^J'^z) = Gd{z) and G^fiz) = G/(z) on each lattice 
site. 

The FK atom can be solved by using the interaction 
representation, such that the time dependence of opera- 
tors is defined by the atomic Hamiltonian, 

H,t = -^^dU+{Ef-^l)pf + udUpf, (3) 

and the time-dependence of the state vectors is governed 
by an evolution operator which is defined by the A-field. 
For purely imaginary times, f € (0, — the evolution 
operator is 

^(f. A) = Tre+fodrfo^df"Hf',f")dHf')d{f")^ 

where Tf[- ■ ■] orders all the operators within the bracket 
with respect to the position on the line (0, — of their 
time-argument in such a way that the time-arguments 
which are closer to zero appear further to the right and 
there is an overall change of sign whenever the time or- 
dering interchanges two fermionic operators. The time 
evolution of the operators between an initial time f' and 
the final time f is determined by Hat as, 

0(f) = e^(^~^')^»'0(f')e-'(^-^')"'". (5) 

The integration is along the imaginary-time axis (f and f 
are purely imaginary), i.e., we have made Wick's rotation 



f = —ir with respect to the real variable used by Brandt 
and Mielsch^. 

The atomic Hamiltonian Tiat and the time-dependent 
field A(f, f ') define the partition function of the FK atom, 

ZatiX,^i,|3) = Tr<j/T^[exp{-/?Hat}^(/3, A)], (6) 

where the trace is taken over the atomic d- and f-states. 
Since the number of f-electrons is a conserved quantity, 
we can write 

ZatiX) - Zo(A, m) + e-'3(^/-^)Zo(A, ^l~U). (7) 

with Zo(A, /i) the partition function of a d-electron sub- 
ject to the A-field in the absence of f-electrons (n/ = 0). 
That is, Zq{\,^) = TTdTf[exp{-(3no}S{P,\)], where 
Ho = —fjid^d. 

The field A(f,-f') gives rise to fluctuations in the 
d-occupancy, which correspond in the equivalent lat- 
tice problem to the local d-fluctuations due to the d- 
electron hopping. We choose A('f, f') such that it satis- 
fies the same (anti-periodic) boundary condition as the 
imaginary-time Green's function (see below) and expand 
it in a Fourier series along the imaginary axis (recall f is 
purely imaginary), 

A(f-f')-T^e-"(^-^')A„. (8) 

n 

Here, ujn — ■7T{2n + 1)T is the Fermionic Matsubara fre- 
quency and we set ks = 1. The A-field in the com- 
plex w-plane, which is needed for the f-electron propaga- 
tor, can be determined by an iterative procedure using 
the DMFT self-consistency condition for the d-electron's 
Green's function. 

The d-electron Green's function is defined as a func- 
tional derivative of the atomic partition functioni, Gd{f— 
f') — (51nZat/(5A(f', f), which gives, 

Gd(f - f') = -^TidfT, [e-'^«-5(/3, A)d(f)dt(f')] . 

(9) 

G/(t — r') is periodic on the imaginary time axis with 
period 2i/3, and is antiperiodic mod {i(3). It depends 
on the difference of the time-arguments, because we are 
in thermal equilibrium and the system is time-translation 
invariant; it has a discontinuity at f = f ' and is therefore 
a non- analytic function of f — f '. 

To find Gf we define the effective- medium Green's 
function, G°(f — f') = (51nZo/(5A(f',f), which reads 

Gdir-f) = -^-^Tr,T, [e-^^o Si(3, \)dif)dHf')] , 

(10) 

and satisfies the equation of motion [EOM)i 

y df[(-a^+iM)(5e(f-f)+A(f-f)]G2(f-f') - <5e(f-f') 

(11) 
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with 5c{f) the delta function defined on the fine seg- 
ment (0, — i/3) with the normahzation J dfdc{f) = 1 (i.e., 
Sc{f) — iT^„ exp[u;„f]). Since also satisfies the 
usual periodic boundary conditions, Fourier transform- 
ing the EOM gives the solution. 



and the full Green's function follows as^ 

Wq , Wi 



Gd{iuj„) 



(12) 



(13) 



with Wq and wi the statistical weighting factors of the 
unoccupied and the occupied f-states, respectively. These 
are given bj*^ wi = 1 — wq and wq — Zq/ Zat , where 



Zo(A,A.) = 2e''-/^n 



(14) 



Defining the self-energy function of the FK atom S^; by 
Dyson's equation. 



^d. = [G"] ^ - [Gd 



(15) 



which holds on the imaginary and the real frequency axes, 
we can find the A-field and the statistical weighting fac- 
tors Wo and wi by an iterative procedure. We start with 
a trial self-energy on the imaginary axis and calculate Gd 
from Eq. get [G°]~^ from Eq. lfT31) . calculate wq and 
wi using Eqs. H14I) and 10, recalculate Gd from H13() . and 
find the new from Eq. (|15|l . Once the procedure is con- 
verged on the imaginary axis, we have the weights Wq and 
wi, the chemical potential fi, and the self-consistent solu- 
tion for the A-field. Since wq and wi are just numbers, we 
can analytically continue Gd{iuJn) given by Eq. H13|l from 
the imaginary axis in the complex frequency-plane and 
repeat the iterative procedure to find retarded quantities 
and the spectral function of the A-field. The knowledge 
of the A-field everywhere in the complex w-plane is a nec- 
essary input to find the Green's function of the localized 
f-electrons. 

The imaginary time Green's function of the f-electrons 
is defined by an expression analogous to Eq. (PJ 



Gf if 



1 

Zat 



-Trd/T, [e-''«-^(/3,A)/(f)/t(f')], 

(16) 

it satisfies the same boundary condition, and has the 
same analytic structure as Gd{f — f'). As usual, we ex- 
press G/(f — f') in terms of two analytic functions of 
f — f ' , such that 



G/(f- 



G>(f-f') forlm(f-f') <0 
G<(f-f') forlm(f-f') >0 



(17) 



To find Gj and G^ we introduce the real variables t = it 
and t' — if' , and represent Gf{f — f') by the Matsubara 
sum. 



Gf{f-f')=T 



>GfiiLUn) (18) 



For T — t' G (— /3, 0), which corresponds to Im(f — f') G 
(0,(3), we replace the sum over Matsubara frequencies by 
a contour integral in the complex frequency-plane. Using 
Res/(z)|^^ = — 1//3, where /{oj) is the Fermi function, 
we obtain 



/oc 
du;f{tj)e 
-OO 



A{cu), (19) 



where A{uj) = — ImG^ ((jj)/7r is the spectral function. 
Here, we used the fact that Gf{z) coincides with the 
retarded Green's function, Gj{z), in the upper half- 
plane and with the advanced Green's function, G^{z), 
in the lower half-plane, and that the discontinuity of 
Gf{z) across the cut along the real w-axis is given by 
[Gf{Lj + iO+)-Gf{Lj-iO+)]/2Tri = -ImGf(w)/7r. Since 
T — t' G (— /3,0), the integrand is well behaved for 
Lu — > ±oo (the cutoff at uj — > oo is provided by the 
Fermi function) and the integral defines an analytic func- 
tion (it has derivatives to all orders). Reinstating the 
imaginary times f and f', we obtain 



G<(f-f') = 



dujA{uj)f{u;)e 



-iuj{T — T') 



(20) 



which can be used to perform the analytic continuation 
from the imaginary to real times, f — > t, f' — > t' . 
Defining the Fourier transform of the real-time Green's 
function as, 



G<{t 



duje^M'~''^G<{u;) 



it follows that 



G<{u;) = 2TTA{u;)f{Lj); 



(21) 



(22) 



hence our definition for G^ and G^ is missing a factor 
of i from the standard definition, but agrees with that of 
Brandt and Urbanels^. Similarly, for t — t' G (0, (3) and 
Im(f - f') e (-/9, 0), we use Res/(-z)|^^^ =1//? and 
express the Matsubara sum for G^(t — t') as, 

/oo 
dojA{uj)f{~uj)e-'^(-^-^'\ (23) 
-oo 

This integrand is also well behaved and shows that Gj 
is an analytic function of imaginary times. The analytic 
continuation to the real axis, f — > t, f' — > t' , gives 



/oo 
dujA{uj)[f{Lu) 
-oo 



j(t-t') 



so that the Fourier transform of G^{t) reads, 



G>icj)^27rA{cj)[f{cj) 
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(24) 



(25) 



At half-filling, where A{ui) = A{-uj), we use G^{t-t') = 
[Gy (i' - t)]* and obtain from the inverse of Eq. (j2i|) the 
result. 



2 

A(w) = / dtcoa{u}t)YleG>At). 

Jo 



(26) 
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Thus, the time-ordered Green's function at real times can 
be written as, 

Gf{t-t')^ [ dLjA(cj)[/(cj)-e(t-t')]e~'"^*~*'\ (27) 



0) 



1 and 



with Q(x) the unit step function Q(x > 
e{x < 0) = 0. 

However, these formal manipulations, which reveal the 
analytic properties of the Green's function, do not explic- 
itly provide the spectral function, A{lo), which is needed 
to find Gy or Gj. On the other hand, the imaginary- 
time formalism provides numerical results for the Green's 
function at the Matsubara frequencies^iii but it does 
not reveal the full analytic structure and it does not pro- 
vide the spectral function. Thus, the real-time Green's 
function of the f-electron cannot be inferred directly from 
these formal and numerical results but has to be calcu- 
lated separately. 

To obtain the real-time properties of the f-electron we 
define the contour-ordered Green's function in the inter- 
action representation as 



1 



where 



(28) 



is the analytic continuation of the evolution operator in 
Eq. I0J from imaginary times to times on the contour 
which is depicted in Fig. ^ for the case t' > t. For t < t' 
the contour starts at t runs to t' , goes back to t, and 
ends at t — i(3. Once again, this notation is missing a 
factor of i in the exponent of the evolution operator and 
as a prefactor for the Green's function, from that used 
in standard approaches, but it agrees with Brandt and 
Urbanel^. 



t 



J t 



t -iB 

FIG. 1: Kadanoff-Baym contour for evaluating the equilib- 
rium Green's function for t < t' . For t > t' the contour starts 
at t' and runs to t, then goes back to t' and ends at t' — i/3. 



The time dependence of the operators on the contour 
is defined as (assuming F precedes t on C), 



0(f) 



(30) 



The Tc [■ ■ ■ ] orders all operators within the bracket with 
respect to the position on the contour of their time argu- 



ments, such that, 

T [f m fUP\] - / /(^/^(^') ^' preceding t on C 
[J[i}J IE jj - -f\i')J{t) i preceding t' on C ' 

(31) 

and similarly for the d operators. 

The A-field on the contour is obtained by an analytic 
continuation from the imaginary time axis. Since the 
A-field, considered as a function of imaginary times, sat- 
isfies the same boundary condition as the Green's func- 
tion, it is composed of two analytic pieces which can be 
continued from imaginary times to complex times on the 
contour C . Starting from the Fourier transform given by 
Eq. (IHJ we obtain, in analogy with Eqs. I^njl and if^ . 
the results, 



A<(f-f') = -- 



du;/(w)e-"^(^-^')lmA-^(w). (32) 



and 



1 r°° 

A>(f - f') = / Mf{^) - l]e-'^''^-^'huY\^{Lj). 

I" J-oo 

(33) 

where ImA^(w) is the spectral function of the A-field. 
These integrals define analytic functions of imaginary 
times and can be continued to the contour shown in Fig.Q] 
simply by substituting f — > f and f' — > ?, where i and t' 
are on the contour C. Thus, the contour-ordered A-field 
can be written as, 



d(jjImA^(w) exp[— ia;(t — t')] 



(34) 



where Oc = if f precedes f on C and Oc = 1 otherwise, 
and t and t' lie anywhere on the contour. Restricting t 
and t' to the vertical part of the contour, and using the 
anti-periodicity along the imaginary axis, we can make 
the inverse Fourier transform 



I I dfAc(f,0)e-'^"^ 
/o 

I I c/fAc(f,-i/3)e"'^"^ (35) 
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and recover the spectral formula 



\{iuJn) 



duj 



Im A^(w) 



(36) 



Contrary to the arguments of Gf{t — t'), which define 
the boundaries of the horizontal piece of the contour, 
the arguments t and t' of the contour-ordered A-field can 
be anywhere on the contour C, so that the dynamical 
mean field connects the real and the imaginary parts of 
the contour. In DMFT, the A-field originates from the 
electron hopping on the lattice and it is responsible not 
only for the propagation of particles in real times but for 
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the thermalization of the system as welL Contrary to 
most non-equihbrium problems, in which the real field is 
switched on at some time to and the integration along 
the vertical part of the Kadanoff-Baym contour can be 
neglected as one approaches the steady state, in DMFT 
problems it is essential to integrate over the whole con- 
tour, because the hopping on the lattice (which gives rise 
to the A-field) is always present. 

We can now find the contour-ordered Green's function 
using the Kadanoff-Baym EOM methods, and we con- 
sider the case t > t' . The fermionic operators on the 
contour satisfy the EOMs, 



z--dt(f) ^ _ Up{t)f{t)]d\t), (37) 



and 



tj-f{t) = [{Ef + Ud\t)d{t)]nt)- (38) 

Since f'^{t)f{t) — /^(0)/(0) does not change with time 
(i.e., it commutes with the Hamiltonian), Eq. ()37|l has a 
simple solution, 

d\t) = e'(A'-t^)(f-*')/t(t')/(t')dt(i')^ (39) 

where t' is the initial time on the contour. The evolution 
of f-electrons is more complicated, because the number 
of d-electrons fluctuates in time and we can only write 
the solution as a contour-ordered product. 



/(?) = e-'(^/-^)(*-*')5^(xri/(t'), 



(40) 



where 



S'Axt) = Te exp I j dtj dt'xiil t')S (l)d(F) j , (41) 
and 

X#, F) = -iUQS- ^SS- I'). (42) 
Thus, we obtain 

G}{t-t') = 

X Tr,; [e-'^«-5e(A,)/(t')/nt')] , (43) 

where Sc{Xc) is the evolution operator in the presence of 
the modified time-dependent potential which is due to 
the fluctuation in the number of d-electrons during the 
propagation of an f-electron from the initial time t' to 
the final time t, 



(44) 



Note, the operator sequence f{t')f'^{t') commutes with 
S'c(Ac) and removes all the occupied f-states from the 



trace, so that the f-propagator can be expressed in terms 
of a partition function of an effective d-electron, 



G>(t-t') = 



-2(Ai,Ac), (45) 



where Tio defines the dynamics of a d-electron when there 
are no f-electrons and 



Z{K) = Tr^ \e~P^-S,{K) 



(46) 



is the partition function of such an electron subject to the 
effective Ac-ficld. Because of time-translation invariance, 
we set t' — Q from now on. 

To find the partition function of a d-electron driven 
by the time-dependent Ac-field we use again functional 
derivative techniques to define an auxiliary Green's func- 
tion. 



1 SZ{K) 
Z(A) 5\{t',ty 



(47) 



such that. 



-^Tr, 



Z(A) 



T,e-P"<^d{t)d\t')S,{K)\ . (48) 



Note, the operator dynamics on the contour is now de- 
fined by Ho, 



0(i) =e*(*-*')«oO(^')e■ 



-^(t-t')•Ho 



(49) 



Next, we introduce an auxiliary contour-ordered Green's 
function for a d-electron driven by the x*"field in the 
absence of the Ar-field, 



5S(*"?) 



5\nZo{xt 



where 

Zoixt) = Trrf [T^e-^^^S^ixt)] = 1 



„|3^l-iut 



(50) 



(51) 



is the effective partition function of such a system. Func- 
tional differentiation gives. 



1 



Mxt) 



Tr. 



Tce-^''od4t)d\i')Sc{xt) 



(52) 

The evaluation of -Zo(Xi) and is straightforward (for 
details see R.efs.lTllandlTI). because the evolution operator 
Scixt) does not change the number of d-electrons and the 
Hilbert space for the d-states comprises only two states 
{ud = and Ud = 1). The Green's functions g and 
go depend explicitly on the contour-times t and t', and 
implicitly on the external time t (recall we set t' = 0). 

Taking the time derivatives of g^{t, F) and (?§(t, F) with 
respect to t we find, using Eq. the EOMs, 



dF[gT\UW,t) = S^it- F), 



(53) 
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where 

[gTHlt') = [gor'im + Ac(i" F), (54) 

and 

[g^,r\t, t') = (-|_ + ^^J)5S- t') + Xtil t'). (55) 
In operator form, this can be written as, 

[5']"V = 1, (56) 

where the unit operator 1 has the matrix elements Sdt — 
i'). The Dyson equation for the integral operator g'^ can, 
thus, be written as, 

[5l"' = k1"'[l+5SA.], (57) 

with the operator product implying an integration over 
the contour C . From the definition of g^ in terms of the 
functional derivatives of the partition function Z(\c), it 
follows that, 

Z(A,) = e^"^ [5^r\ (58) 

where the continuous trace of a contour-ordered operator 
is given by the line integral over the contour C. Using 
Dyson's equation, this can be written as, 

Z(Ae) = Zo(xf)eT^ ind+sSA.)^ (59) 

To approximate the continuous trace by a discrete one^i 
we expand the logarithm, 

ln(l + gSA,) = ^i(5SA,)", (60) 

n 

and replace each contour integral by a discrete sum, using 
a discrete quadrature rule 

N 

dti{t) ^Y,wa{u) (61) 

•^'= i=l 

with weights Wi for the discrete set of times {^i, . . . , t^v} 
on the contour C. Then, the multiple integrals in Ea.lGUI 
reduce to matrix multiplication, and we can use the usual 
expression from linear algebra, 

Tr In (1 + .ggAc) = InDet (1 + ^gAJ (62) 

where Det represents an N x N determinant, we obtain 
the final result, 

G>{t) = -ie-'(^/-'^)*Zo(Xt) (63) 

k 

which we calculate numerically (l/A^c is the approxima- 
tion to the delta function on the contour C with Aic the 
width of the interval that includes the delta function; 
for a [midpoint] rectangular quadrature rule, one takes 
W, = Ate). 



III. COMPUTATIONAL RESULTS 

The numerical evaluation of the /-electron Green's 
function appears to be a rather straightforward proce- 
dure: one decides on a step size for the real-time axis 
Apical and for the imaginary-time axis Atjmag of the 
Kadanoff-Baym contour and then calculates out to the 
largest time that is feasible within the limitations of 
the computational resources. In the results presented 
here, we take Aijcai to range from 0.1 to 0.0125. We fix 
A^imag = 0.05. The cutoff in time is always taken to be 
no larger than 80. In order to calculate the /-electron 
Green's function, we need to calculate the determinant 
of a discretized matrix operator. This is done by first 
diagonalizing the matrix, and then taking the product 
of all of the eigenvalues. This step is the most time- 
consuming step of the calculation, because the matrix 
is a general complex matrix, with no special symmetries, 
and the eigenvalues are usually complex- valued (the max- 
imal matrix size that we consider is about 2100 x 2100). 
Since each time t chosen to evaluate the Green's function 
requires a new contour, the grid of points on the time 
axis, where (t) is generated, need not use the same 
spacing as the discretizing grid of each Kadanoff-Baym 
contour used for discretizing the continuous matrix op- 
erator. Usually, we use a time-grid spacing of 0.2 or 0.1, 
because the Green's function does not normally have os- 
cillations that are on a finer scale than that on the time 
axis. Once the Green's function has been calculated on 
the time-axis grid, we perform a Fourier transform to cal- 
culate it on the real frequency axis. We first spline our 
real-time data (using a shape-preserving Akuba spline) 
onto a real-time grid that is twenty times smaller than 
the originally chosen time grid spacing. Next, we nu- 
merically sum the (cosine) Fourier transform of the real 
part of the Green's function to determine the spectral 
function (which is possible only at half filling; for other 
fillings the analysis is more complicated^) . More details 
of the numerics can be found elsewhereiSi. 

The spectral function satisfies a number of important 
properties. Since our calculation is an equilibrium cal- 
culation (even though we are employing a nonequilib- 
rium formalism), the spectral function is nonnegative and 
the integral of the spectral function over all frequency 
is equal to one. Furthermore, the Green's function on 
the real-time axis approaches — 1 as t — s- and has 
an exponentially decaying (and possibly oscillating) be- 
havior at large times. It also increases quadratically in 
t for small times with a curvature that is independent 
of temperature. Unlike the conduction-electron DOS, 
which is independent^^ of temperature, the /-electron 
DOS evolves^ with T. But, because the value at t = 
[and the first and second derivative of ReG^(t)] is the 
same for all temperatures, we find that the deviations of 
the real-time Green's functions (due to changes in the 
temperature) increase at large times. 

When the DOS develops a gap at low temperature, the 
long-time behavior of the Green's function develops sig- 
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nificant oscillations, with an amplitude that can decay to 
zero very slowly. This creates numerical problems, since 
it implies that the cutoff in time needs to be large in or- 
der to be able to accurately determine the DOS. Indeed, 
we will find that this cutoff dependence limits our ability 
to accurately determine the DOS at low temperature. 




Time t [1/t*] 

FIG. 2: Real part of the (greater) /-electron Green's function 
as a function of time for T = 1 and five different values of U 
on the Bethe lattice. 

We begin our discussion with a plot of the real part of 
the (greater) Green's function versus time for T = 1 on 
the Bethe lattice and for 5 different values of U (Fig.^. 
The values of At,.oai are 0.0125 for [/ = 1, 0.025 for L/ = 2 
and U = 5 and 0.05 for {7 = 2 and U = 2.5. Notice how 
the Green's function appears to have just a smooth ex- 
ponential decay for small U, but as U increases, we first 
see the Green's function assume positive values, and then 
we see that it picks up significant oscillations, whose pe- 
riod decreases as U increases. When U is small enough, 
that there is no gap in the DOS, then we find that the 
long-time behavior is exponentially decaying (with oscil- 
lations entering as the critical value of U for the Mott 
transition is approached). In this regime, we can ex- 
trapolate the results for small time out to large time, 
by fitting the Green's function tail with an exponential 
function, and evaluating that function out to long times. 
This allows us to use a smaller time cutoff (and thereby 
a smaller Aiicai), which becomes increasingly important 
at low temperature in order to maintain high quality in 
the data. 

In Fig.O we plot the logarithm of the absolute value of 
the real part of the (greater) Green's function for T = 0.2 
and J7 = 1.5 on the Bethe lattice. Note how the tails of 
the Green's functions show perfect exponential decay at 
large time. The exponent depends strongly on the dis- 
cretization Atreai of the Kadanoff-Baym contour. But 
because of this simple exponential behavior, we can ap- 
pend extrapolated tails to our more accurate calculations 
(with smaller Aireai) and construct good Fourier trans- 
forms. Unfortunately, there does not seem to be any sim- 
ple way to extrapolate the results to the limit Atrcai 
on the time axis. The problem is that the systematic 
error due to a finite Atreai is not arising from a Trot- 
ter break-up which has a simple error term, but rather 
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Time t [1 /t*] 

FIG. 3: Logarithm of the absolute value of the real part of 
the (greater) /-electron Green's function as a function of time 
for T — 0.2 and U — 1.5 on the Bethe lattice. Three different 
values of Atrcai are shown. 

is arising from the discretization of the continuous ma- 
trix operator. When U is large enough that the tails of 
the Green's functions show oscillatory behavior on top 
of the decaying behavior, we have no simple way to ex- 
trapolate the tails out to large time, and the errors of 
the calculations become larger. [One might have wanted 
to extrapolate In 10^(^)1 on the time axis to A< ^ but 
that cannot be done once {t) crosses zero.] 

When we perform the cosine Fourier transform of the 
real part of G-^{t) to get the /-electron spectral func- 
tion, we first spline the real time data (with an exponen- 
tially decaying tail appended if appropriate) onto a grid 
that is 20 times finer using a shape-preserving Akuba 
spline, and then numerically perform the Fourier trans- 
form. Next we try to extrapolate the Fourier transform to 
A^ioai — *■ 0. This extrapolation is performed using an n- 
point Lagrange interpolation formula (pointwise in w) on 
the different DOS generated for the different Aijoai- Such 
a procedure allows for higher order polynomial approx- 
imations to the extrapolation as more data is included. 
We sometimes find that even though we have data for 
a number of different At^cai values, it is most accurate 
to perform a linear extrapolation for the two smallest 
values of Aircai- We call this extrapolation scheme 6- 
extrapolation. 

We judge the accuracy of our calculations in a number 
of different ways. The first thing we do is to compare the 
moments of the DOS [with or without an extra factor of 
f{u!)] to exact results for those moments (see below). We 
also compare the Matsubara frequency Green's functions 
(generated from an independent program that works di- 
rectly on the imaginary axisLSiii) with the Matsubara 
frequency Green's function generated by integrating the 
spectral formula (with the given DOS) for each Matsub- 
ara frequency 

Gf{iujn) = I duiAfiuj)- — . (64) 

One of the important checks of our numerical accuracy 
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comes from a careful comparison of the calculated results 
with a number of different moment sum-rules of the DOS. 
The sum rules can be derived in a straightforward fash- 
ion: the DOS is first expressed as the imaginary part of 
the Fourier transform of the real-time retarded Green's 
function. By introducing complete sets of states, the time 
dependence of the operators can be expressed in terms of 
the many-body energies of the different states (note that 
because this calculation is performed for the impurity, 
one must include the evolution operator of the A field, 
but since it commutes with the /-electron operators, it 
provides no further complications). These can then be in- 
tegrated over time, and when the imaginary part is taken, 
one gets a delta function in frequency, which allows the 
frequency integral to be performed. Finally, any energy 
factors that remain can be replaced by the Hamiltonian, 
and the sums over the complete sets of states can be per- 
formed. In the end, we are left with operator averages to 
evaluate. These results are summarized in Tabled Note 
that when we perform actual calculations, we always add 
a small shift to the DOS in order to satisfy the zero mo- 
ment sum rule to at least one part in 10^ (typically this 
means adding a shift no larger than 0.005 to the spectral 
function). 

We find that sometimes the (5-extrapolation scheme 
does not further improve the accuracy of the spectral 
function. In that case, it is often more accurate to 
use the result generated with the smallest Atroai- In 
other cases, we find that the exact result for the low- 
est Matsubara frequency Green's function is bracketed 
by the calculation with the smallest Aiicai and the S- 
extrapolation result. In that case, we can average those 
two spectral functions in order to produce better agree- 
ment for Gf{iu}o). We call this extrapolation procedure 
Matsubara-extrapolation. It sometimes can improve the 
accuracy of the results. 

As a general rule of thumb, if we can achieve accu- 
racy of better than 1% for all of the spectral moments, 
and we can achieve four digits of accuracy for all of the 
Matsubara frequency Green's functions, then the result- 
ing DOS is numerically quite accurate. The deviations 
from the exact result are most likely occurring at small 
frequencies, where we need long-time data to get an ac- 
curate Fourier transform, and at high frequencies, where 
the tails don't always decay exactly to zero. 

We illustrate these extrapolation procedures with the 
case J7 = 1.5 on the Bethe lattice (results for [/ = 1 have 
also appearediS). We first focus on high temperature, 
with r = 5. The results for the moment sum rules, for 
the shift to the DOS, and for the lowest Matsubara fre- 
quency Green's function are presented in Table |n] We 
plot the DOS for different Atroai in Fig. 01 As can be 
seen in the figure, as the discretization size decreases, 
the DOS approaches a limiting result, which is close to 
the one predicted by the ^-extrapolation procedure. An 
examination of the table shows how the moment sum 
rules and the Matsubara frequency Green's functions are 
all improved as the discretization error is reduced. The 



extrapolation formula used a quadratic Lagrange inter- 
polation with all the three DOS calculated at different 
Aiicai's. These results show that a systematic extrapola- 
tion procedure is sometimes possible, and that the overall 
accuracy that can be achieved is quite high (of course it 
is difficult to estimate the pointwise accuracy of the DOS 
from any of these integral sum rules). 




-2 -1.5 -1 -0.5 0.5 1 1.5 2 
Frequency w [t*] 

FIG. 4: /-electron DOS for different discretization sizes. 
Also plotted is the 5-extrapolated result using a three-point 
(quadratic) Lagrange interpolation formula. The parameters 
are U = 1.5 on the Bethe lattice, with T = 5.0 

As the temperature is lowered, the spectral function 
sharpens when U is small. In the noninteracting case, the 
spectral function is a delta function. In the interacting 
case, the spectral function approaches a delta function, 
but always maintains a finite width, even at T = 0. Nev- 
ertheless, the calculations become more difficult at lower 
T, because a narrow peak in the DOS implies a slow ex- 
ponential decay in (t), and we find that the discretiza- 
tion error also grows as T is reduced. To illustrate this 
phenomena, we show results for U — 1.5 and T = 0.1. 
The moment sum rules are summarized in Table Hill and 
the DOS are summarized in Fig. |S1 One can see that as 
Aireai IS made smaller, the peak in the DOS is reduced 
in height and increases in width. Furthermore, the 5- 
extrapolation scheme seems to overcorrect, by producing 
a DOS that is too wide (we use a two-point [linear] in- 
terpolation formula here) . The Matsubara-extrapolation 
procedure is much better, but the overall accuracy is re- 
duced relative to the higher temperature results (we find 
only about 3% accuracy for the moments, and three parts 
in 10'^ accuracy for the Matsubara frequency Green's 
functions). We find that this behavior is generic for our 
calculations — usually the calculations are more difficult 
at lower temperature, often requiring a smaller discretiza- 
tion size for the same level of accuracy. We also find 
that the real-frequency extrapolation procedures start to 
break down as T is reduced too. 

A summary of the results, for the case oiU — 1.5 on the 
Bethe lattice is shown in Fig. We have used the most 
accurate DOS calculated at each temperature, by one 
of the two extrapolation procedures. We also included 
the conduction electron DOS, which has a dip develop 
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TABLE I; Sum rules for the /-electron DOS. The expectation value (O) denotes Tr[e"'^'"'"S(A)C']/.Zat. The column on the 
far right gives the results for the half-filled case considered here. The symbol Xmixed denotes the mixed static local charge 
susceptibility between the conduction and the localized electrons. Recall at half filling Ef = —U/2, wi = 1/2, and pe = 1/2. 
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Half-filling result 
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TABLE IL Table of the accuracy of the different calculations of the DOS by comparing results for the different sum rules 


(U = 1.5 on the Bethe lattice, with T = 5). The frequency cutoff for the zero moment sum rule is Ico 


< 15, while all other 


moments are cutoff at the point where the integral stops 


ncreasing and approaches a constant (there is 


usually a decrease for 


larger values of u 


) which normally corresponds to \uj\ < 4 






Moment 


At.cal = 0.1 At.oal = 0.05 


Atrcai ~ 0.025 5-extrapolation 


exact 


1 


0.999994217 1.00000354 


1.0000113 1.00000099 


1 


/(^) 


0.502728929 0.501524381 


0.500528293 0.49994769 


0.5 




-0.023149398 -0.025185717 


-0.026169934 -0.02785382 


-0.027912 




0.465833709 0.507191087 


0.527187957 0.56173037 


0.5625 




0.232916854 0.253595543 


0.263593979 0.28086519 


0.28125 




-0.063584452 -0.063571673 


-0.063545431 -0.06351812 


-0.063518334 


shift 


0.00485 0.005015 


0.002534 -0.0000039 






should increase to about 4.5 as T ^ 0. Note the major 
differences between the localized electron DOS and the 
conduction electron DOS. The /-electron DOS sharpens 
and concentrates much weight around lo = 0, while the 
conduction electron DOS has a dip there. Notice further 
that there is no significant change that we can see in 
our data near the band edge of the conduction electron 
DOS that is also seen in the /-electron DOS, although 
we expect that at T = 0, the bandwidths of both DOS 
should agree with each other. 

We next investigate the case U = 2.5 on the Bethe 
lattice (results for the critical interaction strength U = 2 
appear elsewhere^). This case corresponds to lying just 
on the insulating side of the metal-insulator transition 
(which occurs at [/ = 2). The summary plot of the DOS 
is presented in Fig. [7| Note how the localized electron 
DOS sharpens and develops a gap as T is lowered. What 
is interesting, is that the DOS seems to pile up near the 
correlation-induced gap at low T. We also see a kink 
start to develop near the upper and lower conduction 
band edges, indicating that the /-electron DOS will likely 
vanish outside of the band as T — > 0. Numerically, these 
calculations are challenging. If the discretization error is 
too large, or the time-domain cutoff is too small, then 
we can find negative DOS in the gap region at low tem- 
perature. In fact, the poor quality of our data for larger 
Aireai IS the rcasou why we cannot extrapolate the low 
temperature data faithfully. The accuracy of our calcu- 
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FIG. 5: /-electron DOS for different discretization sizes. Also 
plotted is the 5-extrapolated result using a linear interpolation 
formula for the smallest two Atrcai values and the Matsubara 
extrapolated result. The parameters are U = 1.5 on the Bethe 
lattice, with T = 0.1. 



at the Fermi energy. One can see that the /-electron 
DOS grows and sharpens as T is reduced. We find that 
calculations at much lower temperatures than presented 
here become problematic due to discretization and time- 
domain cutoff errors. Inset into the figure is a plot of 
1/Af{uj = 0) versus T. We have linearly extrapolated the 
last few points to estimate how big the DOS would grow 
as T ^ 0. Our estimate shows that the peak in the DOS 
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TABLE III: Table of the accuracy of the different calculations of the DOS by comparing results for the different sum rules 
{U — 1.5 on the Bethe lattice, with T — 0.1). The frequency cutoff for the zero moment sum rule is \lu\ < 15, while all other 
moments are cutoff at the point where the integral stops increasing and approaches a constant (there is usually a decrease for 
larger values of u) which normally corresponds to < 4. 



Moment 
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At.eal = 0.025 


(5-extrapolation 
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exact 
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1.00000776 


1.00000629 
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1.0000015 
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0.503123093 
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0.5003576 


0.5 
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0.28125 


GfiiuJo) 


-2.08188686 


-1.96833103 


-1.91591427 


-1.86349693 


-1.8890850 


-1.88908508 
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-0.0000091 
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FIG. 6: /-electron DOS for different temperatures. Also 
plotted is the conduction-electron DOS (which is tempera- 
ture independent). The parameters are U = 1.5 on the Bethe 
lattice. Inset is a plot of 1/Af{ijj = 0) versus T. Note how it 
appears to behave linearly at small T allowing us to extrapo- 
late to the T = Q result, so we can predict the maximal height 
of the /-electron DOS at T = 0. 

lations is usually better than 1.5% for the first moment, 
better than 4% for the second moment, and better than 
0.5% for the Matsubara frequency Green's functions (at 
higher temperature, we do significantly better). 

Next we consider a large-gap insulator on the Bethe 
lattice, with U = 5. The summary plot is shown in Fig.|Sl 
Note how the conduction-electron DOS has a large gap 
of about 2.5t*. At high temperature, the /-electron DOS 
has significant subgap states. As T is lowered, we find a 
transfer of spectral weight out of the gap region, with the 
peaks moving towards the gap, and then some additional 
weight being transfered to shoulders that lie close to the 
conduction band edges. The small oscillations in the gap 
region for T = 0.8 are artifacts of the cutoff in time. 

We also present results for the hypercubic lattice. The 



FIG. 7: /-electron DOS for different temperatures with 
U = 2.5 on the Bethe lattice. Also plotted is the conduction- 
electron DOS (which is temperature independent). The T = 5 
data uses 5-extrapolation, the T = 1 and 0.6 data use the 
Matsubara-extrapolation procedure, and the lower tempera- 
tures are not extrapolated (but have At^cai = 0.05). Note how 
the /-electron DOS develops a gap as T is lowered. Note fur- 
ther that a kink starts to develop near the upper and lower 
band edges of the conduction DOS as expected too. Our 
computational accuracy is worst for the subgap DOS at low 
temperature. 



hypercubic lattice does not develop gaps at the metal- 
insulator transition due to the infinite exponential tails 
of the noninteracting Gaussian DOS. But the spectral 
function is suppressed to zero at the Fermi energy and 
there is a "gap region" where the DOS remains exponen- 
tially small. The transition occurs at J7 = ^/2, and we 
expect results for the hypercubic lattice to be similar to 
those of the Bethe lattice when C/Bothc = V^C^hyporcubic 
Brandt and Urbanek's original work* presented results 
for the hypercubic lattice. Unfortunately they gave no 
details on the step sizes used in their computations or of 
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FIG. 8: /-electron DOS for different temperatures with 
(7 = 5 on the Bethe lattice. Also plotted is the conduction- 
electron DOS. No extrapolations are used with this data, just 
the lowest Atroai that produces a positive DOS everywhere, 
and does not develop significant oscillations due to the time- 
domain cutoff. Our computational resources don't allow us 
to go below T = 0.8 here. 



the accuracy of their results. The one discussion of mo- 
ments that they include gives an improper value to the 
second moment of the /-electron DOS, and it is likely 
they never checked the numerical accuracy of their re- 
sults against any moment sum rules. 

We calculate three different values of U for the hy- 
percubic lattice: U = I which has a dip in the conduc- 
tion electron DOS (similar to U ~ 1.5 for the Bethe lat- 
tice), U — 2 which is a "small-gap" insulator (similar to 
U — 2.5 on the Bethe lattice) and U — 4 a, "large-gap" 
insulator (similar to J7 = 5 on the Bethe lattice); the 
near-critical point U — 1.5 appears elsewhere^. Brandt 
and Urbanek showed two DOS for [/ = 1 and five DOS 
for f/ = 2. They did not calculate the U = A case. 

The U — 1 case is plotted in Fig. El The results 
shown here are quite similar to those on the Bethe lat- 
tice (Fig. O. The DOS sharpens as T is lowered, even 
though the conduction electron DOS has a dip at the 
Fermi energy. Inset is a plot of the inverse of the DOS 
at the chemical potential versus T. We can use it to ex- 
trapolate to T = and predict that the spectral function 
peaks at about 4.5. Our results at high temperature and 
at low temperature agree well with those of Brandt and 
Urbanek^. 

Next we consider the case C/ = 2 in Fig. ^| which 
should be compared to the similar results on the Bethe 
lattice (Fig.[7|). Here we see the same kind of behavior — 
the gap is filled at high temperature; as T is lowered, 
spectral weight transfers from the gap region out to the 
band edges; and the peaks of the DOS migrate toward 
the gap regions. Note that the data shown for T — 0.2 
actually has a small region of frequency where the DOS 
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FIG. 9: /-electron DOS for different temperatures with U = 
1 on the hypercubic lattice. Also plotted is the conduction- 
electron DOS. All of the data was extrapolated with one of 
the two extrapolation techniques discussed in the text. Note 
the similarity with Fig. |S|for the Bethe lattice. Inset we plot 
the inverse of the DOS at the chemical potential. Here the 
low-temperature results don't appear to behave in quite the 
linear fashion we saw on the Bethe lattice, but we can still 
attempt to extrapolate to T = with the prediction that the 
peak in the DOS will also be around 4.5 at T = 0. 



goes negative. This is an artifact of the discretization 
error and the time-domain cutoff. 

The results for [/ = 4 on the hypercubic lattice are 
presented in Fig. The behavior is what one expects: 
at high temperature, the gap region is filled in by ther- 
mal excitations. As the temperature is lowered, the gap 
region develops, with spectral weight being transferred 
from the gap out to higher energy. As the temperature 
becomes even lower, the computational needs exceed our 
resources. Note how the peaks in the /-electron DOS are 
pushed closer to the gap region, than the peaks in the 
conduction-electron DOS. 



IV. CONCLUSIONS 

In this contribution we calculated the /-electron DOS 
of the Falicov-Kimball model at half filling. The proce- 
dure requires us to generate the greater Green's function 
for real time and the Fourier transform to get the DOS. 
Unlike the conduction DOS, which is temperature inde- 
pendent, the /-electron DOS has significant temperature 
dependence. For small [/, the DOS sharpens as T is 
lowered to a single-peak structure with a narrow width. 
For large J7, the DOS develops a gap at low T and the 
peaks of the DOS push close to the correlation-induced 
gap edges. When we compare results for similar U val- 
ues on the Bethe and hypercubic lattices, we see similar 
behavior in the DOS. 
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FIG. 10: /-electron DOS for different temperatures with U = 
2 on the hypercubic lattice. Also plotted is the conduction- 
electron DOS. The data is either extrapolated with one of 
the two extrapolation techniques discussed in the text, or we 
work with a fixed value of the discretization on the real-time 
axis. Note the similarity with Fig. |7|for the Bethe lattice. 
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We performed an in depth analysis of the accuracy 
of the numerical calculations. Errors arise from a finite 
discretization error (discretizing the continuous matrix 
operator into a discrete matrix) and a time-domain cut- 
off error [representing the largest time t that G^{t) is 
calculated out to]. We use linear and quadratic moment 
sum rules and the spectral formula for the Matsubara 
frequency Green's functions to gauge the accuracy of the 
calculations. In general, the discretization error becomes 
worse as T ^ and it is quite challenging to get accurate 
results at low temperatures and strong coupling. 

This study is useful to understand problems with the 
accuracy of truly nonequilibrium calculations that use 
similar Kadanoff-Baym contours. While we would not 
have exact sum rules to compare to anymore, it is clear 
that one needs to perform systematic studies in the dis- 
cretization size along the contour to gauge the accuracy 
of the results. One also needs to reduce the real-axis dis- 
cretization as the temperature is reduced. True nonequi- 
librium problems evolve in an external field, and such 
a field can be added into the analysis given here. The 
complicated aspect is being able to construct the local 
Green's function from the local self energy, as the cou- 
pling to a vector potential enters into the hopping part 
of the Hamiltonian, and the local Green's function is no 
longer represented by a simple Hilbert transform. 

We only examined the half-filled case here. This pro- 
vides a significant simplification, as the DOS can be cal- 
culated by a Fourier transform of the real part of (t) . 
For other fillings, the analysis is more complicated and 
usually requires using particle-hole symmetry to gener- 
ate the full DOS. We plan to examine that case in the 
future. 
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FIG. 11: /-electron DOS for different temperatures with 
[/ = 4 on the hypercubic lattice (the T = 5 and T = 2 data 
overlap). Also plotted is the conduction-electron DOS. The 
T = 5 data is calculated with the ^-extrapolation technique; 
all other temperatures work with At^cai ~ 0.05. Note the 
similarity with Fig. |H] for the Bethe lattice. 
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